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ON ! Abstract 

The phenomenon of bursting, in which streaks in turbulent boundary layers oscillate and 
then eject low speed fluid away from the wall, has been studied experimentally, theoretically, 
^L ■ and computationally for more than 50 years because of its importance to the three-dimensional 

'■^ I structure of turbulent boundary layers. We produce five new three-dimensional solutions of 

turbulent plane Couette flow, one of which is periodic while four others are relative periodic. 
Each of these five solutions demonstrates the break-up and re-formation of near-wall coherent 
structures. Four of our solutions are periodic but with drifts in the streamwise direction. More 
surprisingly, two of our solutions are periodic but with drifts in the spanwise direction, a possi- 
bility that does not seem to have been considered in the literature. We argue that a considerable 
^.^1 part of the streakiness observed experimentally in the near- wall region could be due to spanwise 

r^ ■ drifts that accompany the break-up and re-formation of coherent structures. We also compute 

OhI a new periodic solution of plane Couette flow that could be related to transition to turbulence. 

The violent nature of the bursting phenomenon implies the need for good resolution in the 
^\J I computation of periodic and relative periodic solutions within turbulent shear flows. We ad- 

^ ■ dress this computationally demanding requirement with a new algorithm for computing relative 

CN I periodic solutions one of whose features is a combination of two well-known ideas — namely 

the Newton-Krylov iteration and the locally constrained optimal hook step. Each of our six 
solutions is accompanied by an error estimate. 
f~^ . In the concluding discussion, we discuss dynamical principles that suggest that the bursting 

^^3 ' phenomenon, and more generally fluid turbulence, can be understood in terms of periodic and 

^D i relative periodic solutions of the Navier-Stokes equation. 
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1 Introduction 



Turbulent boundary layers are characterized by a viscous sublayer, in which the mean streamwise 
velocity nearly equals the distance from the wall in wall units, a buffer layer, and a logarithmic 
boundary layer. The impressive agreement of theoretical predictions of the me an streamwise veloq 



ity in the viscous sublayer and the loga rithmic boundary layer with experiment (JMonin and Yaglom 



197ll . p. 273) and with computation ( Kim et all 119871 ) can be considered an outstanding success 
in the effort to understand turbulence in fluid flows. It was initially believed that the flow in the 
viscous sublayer was laminar, but experiments in the 50s and 60s lead to the conclusion that ran- 
dom fluctuation s exist in this layer even t hough the mean streamwise velocity in this layer had a 



laminar profile (JMonin and Yagloml . Il97ll . p. 270). 

The phenomenon o f bursting becanie evident during i nvest igations of the viscous sublayer 
and the buffer region ( Klebanoff et al.l . Il962l : iKline et al.l . 119671 ). In this phenomenon, streaks 
in the near-wall region break-up and re-form in a quite striki ng manner. Muc h of the turbu- 
lent energy production occurs in the buffer and viscous layers ( Kline et al.l . 119671 ) and the three- 
dimensional structure of turbulent boundary layers appears to be intimately related to bursting. 
Because of this connection and because of the striki ng nature of the phenornenon itself, burst 



ing has been the subject of numerous experimental (Acarlar and Smithl . 119871 : iBech et al.l . Il995l : 



Klebanoff et al.. 1962; Kline et al 



ical studies (IHamilton et al.l. Il995l: 



20051 : iKawahara and Kida , 



2001 



Holmes et al. 



1967 1_ Smith and Metzleii. 1983 ) and cornputatiqnal or theoret 



19961: lltano and Tohl . l200ll . l2005l : I.Timenez et ah 



Schoppa and Hussainl. 120021') . Th e re has been some discussio n in 



the literature of exactly what is meant by bursting l^ltano and Tohl . l200ll : I Jimenez et al.l . 120051 ). In 



this paper, bursting will always refer to the break-up and re-formation of coherent structures, such 
as streaks, in turbulent buffer regions. 

Although bursting has been much studied, its dyn amics has proved elusive. A large nuinber of 



mechanisms have been proposed to explain bursting (IHam ilton et al.l . Il995l : lltano and Tohl . 12001 



Klebanoff et al.l . Il962l : iKline et al.l . Il967l : ISchoppa and Hussain . .20021 ) . The word mechanism in this 
context does not refer to new physical principles. There is no doubt at all that the incompressible 
Navier-Stokes equation is adequate to explain the dynamics of bursting, and the physics of bursting 
is the physics that goes into that equation and its boundary conditions. However, as is well 
known, the nature of the solutions of the Navier-Stokes equation in the turbulent regime is poorly 
understood. These mechanisms try to provide a way to approxivfiately understand some solutions 
of the Navier-Stokes equation. 

Although approximate, some mechanisms can be useful for computing exact solutions as shown 
by the con s truct ion of traveling wave and steady solutions of channel flows by Waleffe (.Waleffe . 



; 998l. I2OOII. l2003l'l. The me t hod introduced in that bo dy of work was later adapted to pipe flows 
(JFaisst and Eckhardtl . l2003l : IWedin and Kerswelll . 12004 ). However, the break-up and re- formation 



of coherent structures cannot be studied using traveling waves or steady solutions. 

It is certainly very desirable to find exact solutions of the Navier-Stokes equation that corre- 
spond to bursting. These solutions will provide a solid and reliable route to understanding the 
dynamics of bursting. In this context, it is noteworthy that the self-sustaining pr ocess mechanism 



indeed suggested the existence o f periodic solut ions that corr espond to burst i ng (|Hamilton et al. 



I995I : iKawahara and Kidal . I2OO1I : IWaleffel . Il997l '). We follow IHamilton eraP (ll995|) and conduct 
our computations using plane Couette flow at a Reynolds number {Re) of 400. In plane Couette 
flow two parallel walls move in opposite directions with equal speed and drive the fluid in-between. 
The Reynolds number is based on half the separation between the walls and half the difference 
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Table 1: This table gives data for six periodic and relative periodic solutions of plane Couette 
flow labeled Pi through Pq. The markers in the first column are used to distinguish between these 
solutions in later plots. The shifts Sx and Sz are explained in the text. The T jT^ column gives the 
period, with T^ being the period in wall units. The Xmax column gives the maximum characteristic 
multiplier. The last three columns report the frictional Reynolds number, the resolution of the 
computational grid, and the relative error in the computation. 



between the wall velocities. To render the computational domain finite, we assume the domain to 
be periodic in the streamwise and spanwise directions, with periods equal to 27rA^ and 21: K^- Un- 
less otherwise stated, our computations use K^, = 0.875 and A^ = 0.6 to facilitate comparison with 



earlie r computations and because of particular advantages of this box described in (JHamilton et al. 
I995I ). The walls are assumed to be at y = ±1, with the upper wall moving in the x or streamwise 
direction in the positive sense with speed equal to 1. 

It is known experimentally that the details of bur sting in the near-wall regi on are remarkably 
similar over a very wide range of Reynolds numbers ( Smith and Metzlerl . Il983l ). Thus our use of 
Re = 400 in plane Couette fiow is an accept able choice. Turb ulent spots have been observed in 
plane Couette fiow experiments at Re = 360 ( Bech et al.l . Il995l ). 

Table [U reports data for six periodic or relative periodic solutions that we computed. The final 
velocity field obtained by integrating the initial velocity field over one full period equals the initial 
velocity field for periodic solutions. But for relative periodic solutions the final velocity field equals 
the initial velocity field after shifts equal to Sx and Sz in the streamwise and spanwise directions, 
respectively. These shifts, which are reported in Tabled! are both for Pi and Pg. Thus both 
those solutions are periodic. The equations of plane Couette fiow are unchanged by translations 
in the streamwise and spanwise directions. The existence of relative periodic solutions is possible 
because of those invariances. 

Frictional velocity and frictional length can be obtained using the mean shear at the wall 
(|Monin and Yagloml . Il97ll . p. 265). Those quantities are the basis of wall units. Throughout this 
paper, wherever wall units are used, the mean shear is obtained by averaging at the upper wall 
over one single period of a periodic or relative periodic solution. Following standard practice, the 
use of wall units is signaled by using + as a superscript. The frictional Reynolds number Rer was 
obtained using the mean shear at the upper wall and the distance between the two walls. Table [1] 
shows that the Rer for Pi is significantly lower than it is for the other solutions. This is because 
Pi alone is related to tr a nsitio n to turbulence, while all the others are related to bursting. 



Kawahara and Kidal (|200ll ) found a periodic solution of plane Couette fiow with period equal 



to 85.5 which appears similar to Pi. Howev er, their so l ution satisfies the shift-refiection and shift- 
rotation symmetries of plane Couette flow ( Kawaharal . l2005l ). Although Pi has zero mass flux in 



the streamwise direction like their solution, it is very far from satisfying either symmetry as will be 
shown. Therefore Pi is a new solution. The bursting solutions P2 through Pg are all new, and we 
will presently argue that these are the first computations of bursting periodic or relative periodic 
solutions that demonstrably correspond to solutions of the Navier-Stokes equation. 

All the characteristic multipliers Xmax listed in Table [1] are outside the unit circle, which implies 
instability of the solutions Pi through Pg. Yet we report relative errors for all these solutions. 
A relative error of 10~^ implies that the computed solution matches a solution of the Navier- 
Stokes equation up to at least 7 digits. The manner in which these error estimates were found 
is explained in Section 3. Significantly, these error estimates can be verified using any good DNS 
(direct numerical simulation) code in spite of the instability of the underlying solutions. Such 
quantitative reproducibility is a step forward for turbulence computations. 

The numbers 2L, M, and 2N in Table [1] give the number of grid points in the x, y, and 
z directions. We used a Fourier grid in the x and z directions and a Chebyshev grid in the y 
direction. Good spati al res olution is the key to finding solutions that are not numerical artifacts. 
Hamilton et~aD (jl995! ) and ' Kawahara and Kidal (120011 ) used (2L, M, 2N) = (16, 32, 16) in their 



plane Couette flow computations. The evidence for striking recurrences and the existence of periodic 
solutions offered in these works is significant. However, estimates for spatial discretization error 
described in Section 3 imply that the spatial discretization error with (2L, M, 2N) = (16, 32, 16) 
for bursting solutions is at least 5% and can be twice as much. 

A section 2, we describe a new method for finding periodic and relative periodic solutions. The 
number of degrees of freedom that determine the initial velocity field for P3 is 319790 and the 
number of degrees of freedom for P2, P5, and Pg is 273918. These numbers exceed the number 
of degrees of freedom in any earlier computation of periodic solutions by at least a factor of 20, 
and our method can also compute relative periodic solut i ons. One feature of the method is a 
combination of the Newton-Kr ylov iteration (Kellev^ . l2003l : I Sanchez et al.l . |200J) with the locally 



constrained optimal hook step (jDennis and Schnabel 119961 ) . The locally constrained optimal hook 



i 



step is related to the Levenberg-Marquardt procedure and is a well-established idea in optimization. 
However, its possibilities seem to have been largely overlooked in computations of periodic and 
steady solutions. The combination of Newton-Krylov iterations with the locally constrained optimal 
hook step is simple, but powerful. It is much more effective than the often used damped Newton 
iteration. 

In Section 4, we develop the connection of the computations summarized in Table [T] to the 
dynamics of the bursting phenomenon. There has been some discussion about whether the near- 
wall bursting observed experimentall y is due to the advect ion of coherent objects or to the break-up 



and re- formation of coherent objects (I Jimenez et al.l . l2005l ) . The relative periodic solutions reported 



in Table [U and especially the spanwise drifts of P2 and P4, will be shown to be significant in this 
respect. 

In the concluding Section 5, we discuss our belief that a good route to understanding the 
dynamics of a differential equation is by computing its solutions and recognizing the relationship 
between those solutions. We discuss dynamical principles that suggest that infinitely many periodic 
and relative periodic motions can be found within turbulent flows. About half a century ago, a 
common belief was that linearly stable solutions are observed in nature and in experiment, while the 
unstable ones are not. Although all these infinitely many periodic and relative periodic motions 
are certain to be linearly unstable, their instability is only a manifestation of the instability of 
turbulent flows. In spite of their instability, these solutions are relevant both to natural phenomena 



and experiment as we demonstrate in Section 4 and as we argue in Section 5. 

The idea of understanding phenomena using well-resolved and linearly unstable nonlinear so- 
lutions is beginn ing to take root in recent research on the transition to turbulence in shear flows 
( Kerswell 12005). Linearly u nstable nonlinear traveling waves have been observed in a pipe flow 



experiment ( Hof et al.1 . |2004| ). Although transition to turbulence in pipe and channel flows is still 



an unsolved problem, there c an be little doubt that the nonlinear traveling waves such as the ones 
observed bv iHof et alj ( 2004 ) will be an important part of an eventual solution of the transition 



problem. Periodic and relative periodic solutions are the next step for understanding turbulent 
phenomena such as bursting. 

2 A numerical method for finding relative periodic solutions 

This section describes a numerical method for finding periodic and relative periodic solutions in 
plane Couette flow. Although the description of the numerical method is specific to plane Couette 
flow, it can be easily adapted to other partial differential equations. Our method has three new 
aspects. Firstly, we describe a way to find good initial guesses. Secondly, we show how to set up 
the Newton equations for finding relative periodic solutions. Thirdly, we show how to modify the 
Newton-Krylov procedure to compute the locally constrained optimal hook step. 
The Navier-Stokes equation for incompressible flow takes the form 

du/dt + (u.V)u = -(l//))Vp + (l/i?e)Au (2.1) 

with the incompressibility constraint implying V.u = 0. For plane Couette flow, the boundary 
conditions are u = (±1,0,0) at the walls (which are at y = ±1) in the y or wall-normal direction, 
and periodic in the other two directions with periods equal to 27rAa; and 27rA2. The equation 
cannot be viewed as a dynamical system when written in this form because the velocity field u has 
to satisfy the zero divergence condition and there is no explicit equation for evolving the pressure p 
in time. It is necessary to rewrite the equation as a dynamical system for the purpose of computing 
periodic and relative periodic solutions. 

The pressure term can be completely eliminated by recasting ()2.ip in terms of the wall-normal 
velocity v{x,y,z), the wall-normal vorticity r]{x,y,z), and the mean components u{y) and w{y) of 
the streamwise and spanwise velocities. The boundary conditions become n(±l) = ±1, t()(±l) = 0, 
and v{x, ±1, z) = Vy{x, ±1, z) = r]{x, ±1, z) = 0. The velocity field u can be constructed from u, w, 
V, and r] using V.u = 0. The velocity component v and the vorticity component r/ are discretized 
in the x and z directions using Fourier modes. For example. 



' lIx 

-N<n<N 



/ \ X — ^ -^ / N / ilX iTiZ\ . , 

'v{x,y,z)= 2^ 7;i^„(y)exp(^— + — j, (2.2) 

-L<l<L "' '' 



where the dependence on t is not shown explicitly. The Fourier modes t}/^^ are defined by replacing 
V b y 7? in 



Kim et al.l (|l987l ) proposed a very good numerical method for integrating the Navier-Stokes 
equation using this formulation and we follow their approach. The Navier-Stokes equation (|2.ip 
becomes a dynamical system in this formulation. 

We use a Fourier grid with 1L and 2A^ points in the x and z directions, but we always set the 
modes with I = L or n = N equal to 0. We use M + 1 Chebyshev points in the y direction. After 
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Figure 1: This figure shows the projection of a trajectory of plane Couette flow onto the energy 
dissipation (D) and energy input (/) plane 

spatial discretization, the dynamical system has approximately 8LMN degrees of a freedom. A 
count that takes into account the boundary conditions on u, w, v and ry, the treatment of modes 
with I = L or n = N, and the fact that the mean components of v and rj are always zero gives 
2(M— 1)+(2M— 4)((2A^— 1)(2L— 1) — 1) as the exact number of degrees of freedom. It is important to 
get the number of degrees of freedom exactly right. The spatially discretized system can be thought 
of as a dynamical system of the form X = /(X), where the 2(M-1) + (2M-4)((2A^-1)(2L-1)-1) 
components of X encode u, w, v and rj. All the components of X are real numbers. 

Our code implements the nonlinear terms in advection, rotation, and skew-symmetric forms. 
The viscous terms are treated implicitly and the advection terms are treated explicitly when dis- 
cretizing time. The code employs dealiasing using the 3/2 rule in the x and z directions. The code 
was tested using three-dimensional modes of the Orr-Sommerfeld equation and in various other 
ways. 

2.1 Finding initial guesses 

To find an initial guess for a relative periodic solution, we begin by looking at projections of 
trajectories to the energy dissipation and energy input plane. The rate of energy dissipation per 
unit volume for plane Couette flow is given by 



D 



1 



27rAz 



+ 1 

-1 



2itAx 



Svr^A^A^ JO J-i JO 
and the rate of energy input per unit volume is given by 



Vu\ +\Vv\ + |Vw| dxdydz 



Svr^AxA, 



27rA:c 



2^A, Q^ 

dy 



du 
y=i dy 



dxdz. 



^1 



(2.3) 



(2.4) 



For the laminar solution (u, f,w) = (y,0, 0), both D and / are normalized to evaluate to 1. 

The trajectory that starts at A in Figure [1] appears to come close at B. But that close approach 
is not significant because even distant points in phase space can coincide in a two-dimensional 
projection. What we look for is the resemblance of the shape of the trajectory from B onwards 
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with the shape of the trajectory from A onwards. In Figure [H the trajectory that begins at A 
develops a type of protrusion in the upper right part of the figure, but the trajectory that begins 
at B does not. Therefore, A would not be a good initial guess for a relative periodic solution. If it 
were, we would examine velocity fields on the trajectory near B, and shift them in the streamwise 
and spanwise directions to bring them as close to A as possible. The initial guesses for the period 
and the shifts would be determined following that examination. 

2.2 Newton equations for finding relative periodic solutions 

Given the Fourier representation (j2.2p of v{x, y, z), the Fourier representation of v{x + Sx,y, z + Sz) 
is given by 

/ N v^ ft-lsx\ f'-nsz\^ , . / dx Lnz\ . ^. 

v{x + Sx,y,z + Sz)= 2^ expl-— lexpl— — jt;/,„(y)exp( — + -— 1. (2.5) 



-L<l<L 

-N<n<N 



Ax J "V A, J '•"'"' "VA^ A, J 



Define the operators Ti and T2 by 



.-I- / N v^ ''^ - / \ f dx inz\ 

Tiv{x,y,z)= 2^ — z;i^„(y)exp( — + -—1 (2.6a) 



~L<l<L 

'N<n<N 



Ax ''"'"' '\Ax A, 
dx inz' 



T2v{x,y,z)= 2^ —vi^n{y)exp\^— + —j. (2.66) 



-~L<l<L ' " ' ^ 

-N<n<N 



The operators Ti and T2 are infinitesimal generators of the group of translations along the stream- 
wise and spanwise directions. The shift in (j2.5p is given by e^^^e'^''^v{x,y,z). To shift a velocity 
field given by u, w, v, and rj by Sx along x and by Sz along z, we need only apply g'^^^e'^^^ to v 
and T]. To shift a velocity field encoded by Xq, we can convert Xq to u, w, v, r], shift v and rj, and 
then convert back. The shift operation on Xq will be denoted by e^^^e^^^Xo. 

The Navier-Stokes equations for plane Couette flow are unchanged by shifts in the x and z 
directions. In terms oi X = f{X), this property becomes /(e'^^'^^e^^'^^x) = e'^'^'^^e'^'^'^'^fiX). 

To find a relative periodic solution of plane Couette flow, we have to find an initial velocity field 
Xq, a period T, and shifts Sx and Sz such that 

e-'-'^^e-''^^X{T;Xo)=Xo. (2.7) 

We first set up the Newton iteration, although the Newton iteration by itself is entirely inadequate. 
Let Xq, Sx, Sz, T be our initial guess for solving ^M) and let Yq = e'''^'^^ e'^-'^-' X{T; Xq). Then 
the relative error in the initial guess is given by 

||yo-^o||/||^oli. (2.8) 

If we assume Xq + 5Xq, Sx + Ssx, Sz + Ssz, T + 5T to be the solution of (|2.8p that is close to the 
initial guess and linearize about the initial guess, we get 

i6sx){-TiYo) + {6sz){-T2Yo) + {ST)f{Yo) + e'^^^^e"^^^^ ^^j,^^°^ SXo = dXo + (Xo - Yo). (2.9) 

dXo 



The number of equations in the hnear system (j2.9p is the same as the dimension of Xq or 6Xq. 
Three more equations are necessary to have as many equations as unknowns. Those equations are 



6X^{T,Xo) -- 


= 


SX^{T2Xo) -- 


= 


6X^f{Xo) -- 


= 0, 



(2.10) 

where 6Xq denotes the transpose of 6Xq . Since the Navier-Stokes equations for plane Couette flow 
are unchanged by shifts in the x and z directions, shifting Xq in the x or z directions will only shift 
Yo in the x and z directions. The first two equations in (I2.10p above require that the correction 
6Xq to Xq must not have components that shift Xq infinitesimally in the x or z directions. If the 
correction 6Xq slightly advances Xq along the flow induced by X = f{X), the corrected velocity 
field will remain on the same orbit. The third equation in (j2.10p requires that the correction 5Xq 
must have no component along /(Xq). The equations (|2.9p and (j2.10p together constitute the 
Newton system. 

It is convenient to rewrite the Newton system as Ma = p, where a = {6Xq;6sx',Ssz',ST) and 
p = (Yq — Xq; 0; 0; 0) are both column vectors. The structure of M follows from (j2.9p and (j2.10p . 

As evident from (12. 9p . the application of M to cr requires the computation of — ^-' "' 5Xq. 
This directional derivative can be computed using differences to about 7 digits of accuracy, which 
is entirely adequate. 

2.3 Finding the locally constrained optimal hook step 

The dimension of the linearized system M can exceed 3 x 10^ for periodic solutions such as P3, and 
it is impractical to form the Newton system explicitly. The Newton step a can be found, however , 
using a Krylov subspace method such as GMRES which is described in ( Trefethen and Baul . 119971 ) 
for example. 

Often pe riodic solutions and steady states are computed within a bifurcation-continuation sce- 



nario as in (jSanchez et al.l . 12004 ). We are not in such a scenario here and the Newton step by 
itself never leads to convergence. The widely used expedient of damping the Newton step is also 
ineffective. 

The locally constrained optimal hook step is based on the idea that given a radius r within 
which we trust the line arization, the best step a is obtained by minimizing \\Ma — p\\ subject to 



the constraint \\a\\ < r ( Dennis and SchnabeJ 119961 ). The radius of the trust region is varied from 



step to step by comparing the actu al reduction in error followi ng the step with the reduction in 
error predicted by the linearization ( Dennis and Schnabeill996l ). 



We show how to compute the locally constrained optimal hook step within a Krylov subspace. 
Let Qd and Qd+i be matrices with d and d + 1 orthonormal columns obtained by applying GMRES 
with the starting vector p. The columns of these matrices are orthonormal bases for Krylov sub- 
spaces of dimensions d and d -|- 1. Before trying to find the locally constrained optimal hook step, 
we make sure that d is large enough to permit a solution of the Newton equation Ma = p with 
small relative residual error given by \\Ma — p\\/\\p\\- Let Hd+i^d be the upper Hessenberg matrix 
that satisfies MQd = Qd+iHd+i,d- We minimize \\Hd+i^dcrd — Qd+iP\\ subject to the constraint 
||o'd|| ^ 1"- The solution of this minimization problem for ad using the singular value decomposition 
of Hd+i^d is feasible because d is typically a number smaller than 30. Once ad is found, the locally 
constrained optimal hookstep is given by Qd^^d- 



Time Step 



Figure 2: The plot above shows the dependence of the global error, obtained by integrating Pi for 
one full period, on the time step for numerical integrators of 2nd and 3rd order. The 2nd order 
integrator used was Crank-Nic olson-Adam-Bashforth and the third order method was the (4, 3, 3) 



method of (JAscher et alJ . Il997l ). 



2.4 Computing traveling wave solutions 

The method for computing relative periodic solutions can be modified to compute traveling wave 
solutions. Only minimal modifications are necessary. Instead of allowing the period T to vary from 
iteration to iteration, the first modification is to fix T at a value that is definitely smaller than the 
period of any periodic or relative periodic solution. The second modification is to drop the last of 
the three equations in (|2.10p to get a square system of equations for 5Xq, 6s x, and 6s z- The speed 
of the traveling wave in the streamwise and spanwise directions will be given by Sx/T and Sz/T, 
respectively. 

The basis of this approach for computing traveling wave solutions is the observation that the 
initial velocity field of a traveling wave solution is a fixed point of the time T map of the flow, if 
the final velocity field is shifted by appropriate amounts Sx and Sz in the streamwise and spanwise 
directions. The shifts depend upon the wave speeds as indicated in the previous paragraph. 

In some cases, syi nmetries of th e velocity field may imply that the traveling wave must actually 
be a steady solution ( Waleffd . l2003l ) . For computing traveling waves in general we solve for Sx and 
Sz to determine the wave speeds. But if it is known in advance that the wave speeds are zero, we 
set Sx = Sz = and drop the first two equations of ()2.10p . We moved some of Waleffe's traveling 
wave solutions ( Waleffd . l2003l ). whose data are posted publicly, from a 32 x 34 x 32 grid to finer 
grids, with (2L, M, 2N) = (48, 73, 48) or better, and then refined the solutions to better accuracy. 
During the refinement the period T was fixed at 1 or 5 or 10, with larger values of T implying 
faster convergence of the GMRES iteration. The refined solutions too were invariant under the 
shift-reflection and shift-rotation symmetries. 



3 Verifiability of computed relative periodic solutions 

Although the solutions reported in Tabled] are all linearly unstable, they can all be verified using 
a good DNS code for channel flow as will be shown in this section. We also explain how the error 
estimates reported in the last column of Tabled] were derived. 
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Figure 3: The first two plots show the variation of the energy with streamwise and spanwise 
wavenumbers, respectively, for the relative periodic solution Pq. The energies are computed using 
the slice y = and at eight equally spaced instants along Pg's period. The third plot shows the 
magnitude of the energy in wall-normal Chebyshev mode against the Chebyshev mode at eight 
equally spaced instants along Pqs period. 



The characteristic multipliers Xmax reported in Table [T] are all greater than 1 but less than 100 
in magnitude. For characteristic multipliers in that range, the errors due to time discretization 
can be made negligible. To see that, consider the initial value problem x = f{x), x{0) = xq. 
The global error after time T is defined as ||x(r;xo) — a;(T;xo)||, where x{T;xo) is the numerical 
approximation to the solution x{T; xq) at time T. For an integrator of order r the global error is 
asymptotically equal io E{T)N' in the lim it of small h. Indeed, explicit formulas can be obtained 
for the function E{T) ( Viswanathl . l200ll ). Those formulas indicate that E{T) will increase with 
Xmax for the solutions listed in Table [TJ However, as indicated by the asymptotic formula and 
as demonstrated in Figure O the global error can be made quite small by taking a small time 
step. For all the six solutions in Table [U we carried out computations similar to the one shown 
in Figure [2] and chose a time step small enough to make the time discretization errors irrelevant. 
The final com putations were all c arried out using a 3rd order implicit-explicit Runge-Kutta method 
developed bv lAscher et al.l (jl997l ). 

The argument in the previous paragraph would have been invalid if some characteristic multi- 
pliers were too large in magnitude. For instance, if Xmax > 10^^, even the rounding errors would 
be amplified to an 0(1) magnitude over a single cycle. In such situations, one has to use multiple 
shooting. 

Thus if the system X = f{X) is obtained by spatially discretizing the Navier-Stokes equation 
(j2.ip and we compute a relative periodic solution of that system with a small enough time step, how 
close the computed solution is to a true periodic solution of the Navier-Stokes equation is entirely 
determined by the spatial discretization error. 

We estimate the spatial discretization error in two ways. The first way is to graph energy against 
streamwise wavenumber, spanwise wavenumber, and wall-normal Chebyshev mode as shown in 
Figure [3l The first two plots in that figure do not show the energy that corresponds to the 
wavenumber 0. Because the Chebyshev polynomials are not orthogonal with respect to the Lebesgue 
measure, the decomposition of the energy into Chebyshev modes is necessarily somewhat arbitrary. 
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We defined 

(■2-kKz r2iT\x 

E{y) = I I u{x,y,zf + v{x,y,zf + w{x,y,zf dxdy, 

Jo Jo 

and expressed E{y) as a linear combination of Chebyshev polynomials to get the third plot in 
Figure El More specifically, if E{y) = J2m=o '^"^'^^(v) ' where Tm{y) are Chebyshev polnomials, 
the third plot in Figure [3] plots \cm\ against m. Estimates for the discretization error are obtained 
by taking the square root of the fraction of the energy in the highest mode. Thus we will have 
streamwise, spanwise, and wall-normal estimates at each point along the relative periodic solution. 
The worst of these estimates corresponds to the direction that is less well-resolved than the other 
two, and the instant along the solution where the velocity field is hardest to resolve. 

Another possibly more reliable way to estimate the spatial discretization error is to take the 
initial state of the computed relative periodic solution and then move it to a much finer grid. The 
initial data is integrated on this finer grid using a sufficiently small time step for one full period. 
The final state is then shifted using the shifts Sx and Sz- The quantity 

1 1 shifted final state — initial state 1 1 



1 1 initial state || 

is taken as an estimate of the spatial discretization error. 

In all six cases reported in Table [TJ these two methods gave comparable estimates for the spatial 
discretization error. All the errors reported in Tabled] were obtained using the second method and 
a finer grid with (2L, M, 2A^) = (64, 90, 64). As our discussion of time discretization error makes it 
clear, all the error estimates can be verified using a good DNS code with a sufficiently small time 
step. Gibson used data abo ut Pi tha.t had a relative error of 10~^ and verified that error estimate 
using his Channelflow code ( Gibsonl . l2002l ). 



4 Relative periodic solutions and the bursting phenomenon 

The most striking thing a bout the bursting phenomenon in experiments is its recurrent nature. As 
Acarlar and SmithI (j 198 71 ) state "The study of boundary-layer turbulence in the last thirty years 



has clearly demonstrated that the chaotic behavior referred to as turbulence has a systematic 
organization. Most researchers share a common belief that a cyclic bursting phenomenon is the 
predominant mode of turbulence production." In this section, we demonstrate the relevance of the 
periodic and relative periodic solutions listed in Table [J to the bursting phenomenon. 

4-1 Position of relative periodic solutions in phase space 

A good preliminary idea of the rel ative location of solut i ons i n phase space can be formed by 



projecting the orbits to the D-I plane (jKawahara and Kidal . 120011 ). The normalization used for D 



and / in (12. 3p and ()2.4p is such that the laminar solution of plane Couette flow is located at (1, 1). 
We immediately see from Figure H] that Pi is much closer to the laminar solution than the other 
solutions. As mentioned earlier. Pi is not a bursting solution, while the others are. 

Greater energy dissipation is connected with steeper gradients in the flow. It is therefore 
tempting to conclude that orbits that travel farther into the upper-right corner of Figure H] are 
harder to resolve. That conclusion is true only approximately. Plots of energy against wavenumber 
or Chebyshev mode look very similar to Figure [3] for P2, P3, P4, and P5, but not for Pi. The 
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Figure 4: Energy dissipation (D) and energy input (/) are defined by (12. Sh and (|2.4p . The periodic 
orbit in the lower left corner is Pi . The other orbits correspond to the solutions P2 through Pg with 
the correspondence given by the markers. The marker for each periodic orbit is listed in Tabled) 
A random turbulent trajectory is shown in the background. 



most striking aspect of the plots for Pq in Figure [3] lies in the plot of energy against streamwise 
wavenumber. While at one instant the energy falls by more than a factor of 10^'^ over the range of 
streamwise wavenumbers represented in the computational grid, it falls by only a factor of about 10^ 
at another instant. This striking change in energy distribution is directly related to the breakdown 
of streaks and is observed in P2 through Pq but not in Pi . 

The existence of periodic and relative periodic solutions is related to t wo general principles in 
dyna mics, namely the Poincare recurrence theorem and the closing lemma (JKatok and Hasselblattl . 
I995I ). The recurrent nature of the bursting process was observed in direct numerical simu- 
lation of plane Couet t e turbu l ence and motivated the derivation of the self-sustaining process 
( Hamilton et al.l . Il995l : IWaleffd . 119971 ). It must be pointed out, however, that the right notion of 
recurrence in plane Couette flow is not that of periodicity but that of relative periodicity. This is 
because if an initial velocity field for plane Couette flow is shifted in either the streamwise or the 
spanwise direction, the final velocity fields after a certain interval of time must be equal to each 
other after exactly the same shifts. From the point of view of the dynamics, velocity fields that are 
related by streamwise and spanwise shifts are equivalent. 

The same general principles that suggest the existence of periodic and relative periodic solutions 
that correspond to bursting, also suggest that there will be infinitely many of them. Indeed, these 
solutions probably give the best or only means to uncover the "systematic organization" found 
within turbulent boundary layers that Acarlar and Smith refer to in the passage quoted at the 
beginning of this section. This possibility will be discussed at length in the concluding discussion 
section. 

4.2 Near-wall statistics 



mic 



Turbulent bound ary layers can be divided i nto a viscous sublayer, a buffer layer, and a logarith- 
boundary layer ( Monin and Yagloml . Il97ll ). In the viscous sublayer, which is approximately 5 
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Figure 5: The first of the top two plots shows the dependence of the mean streamwise velocity 
<u>^ in wall units upon the distance y"*" from the upper wall, also in wall units. The second 
plot shows the dependence of the turbulent intensity as given by the root mean square streamwise 
velocity u'^ upon the distance from the wall. The bottom plot shows the dependence of turbulent 
energy production upon the distance from the wall. The correspondence to the solutions Pi through 
Pq listed in Tabled] is shown using markers. 
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wall units thick, the mean streamwise velocity nearly equals the distance from the wall in wall units. 
The buffer layer extends from 5 wall units to approximately 25 or 30 wall units, and the logarithmic 
boundary layer begins after that. The viscous stresses dominate in the viscous sublayer, while the 
Reynolds stresses dominate in the logarithmic boundary layer. Both stresses are significant in the 
buffer layer. 

The first plot in Figure [5] shows the dependence of the mean streamwise velocity upon the 
distance from the wall. The averages are computed over one single period for each periodic solution 
in all plots in Figure O The law of the wall is reproduced correctly in the viscous sublayer, but 
the logarithmic boundary layer is not fully developed. The distance between the two moving walls 
in plane Couette flow in only about 70 wall units or so for the solutions listed in Table [H which is 
one reason the logarithmic boundary layer is not fully developed. Another related reason is that 
the frictional Reynolds numbers listed in Table [1] are too low for the logarithmic boundary layer 
to be fully formed. A dynamical investigation of the interac tion of structures awa y from the buffer 
region with structures in the buffer region can be found in ( Itano and Tohl . l2005l ). 

The most important features of bursting are in the buffer layer and these features are reproduced 
correctly by P2 through Pg as shown by the other two plots in Figure [5j The top-right plot 
graphs the turbulent intensity against the distance from the wall. The graphs for P2 through 
Pq all have the right shape and the turbulent intensity peaks between 10 and 20 wall units for 
each of those solutions. However, t he peak is slightly e l evate d as compared with the theory and 
experiments recorded in Figure 27 of iMonin and Yaglqrn ( 19711 ) or with the "corrected" experiment 



and computation recorded in Figure 7 of (JKim et al. 
Reynolds numl j er eff ect . 



19871 ). 



(j Jimenez et al.l . l2005l ) . 



This elevation of the peak is a low 
The peak of the turbulence intensity compares well with Figure 3 of 



The energy balance equation, which is obtained by applying the method of averaging to th e 
Navier-Stokes equation, is important both in theory and in practice ( Monin and Yagloml . Il97ll ). 
Physical interpretations can be associated with the terms of that equation, and we will look at the 
so-called turbulent energy production term. This term equals 



<u V > 



d <u> 



dy 



where u* = u— <u> and v* = v— <v> are the fluctuating components of the streamwise and 
wall-normal velocities and <u> is the mean streamwise velocity. In Figure [5l the turbulent energy 
production is expressed in wall units, al though experimen talists do not always use wall units for 



expressing turbulent energy production (JKline et al.l . 119671 ) 



Turbulent energy production has a sharp peak in the buffer region for each of P2 through Pq, 
as shown by the bottom plot in Figure [5j Turbulent energy production can be readily measured 
in ex periments and its sharp peak in the buffer region has intrigued experimentalists for a long 



time (JKline et al.l . 119671 ). The significance of the bursting phenomenon is in part because of its 



connection to turbulent energy production. We observe from the bottom plot of Figure [5] that for 
Pi, which is not a bursting solution, turbulent energy production attains its maximum value farther 
away from the wall. 



4.3 Break-up and advection of coherent structures 

Streaks are the most prominent coherent structures observed in turbulent boundary layers. 
In streaky velocity fields, the streamwise velocity is only weakly dependent on the streamwise 
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Figure 6: Isolines of the streamwise velocity in the z vs.x plane at y^ ~ 10 for the relative periodic 
solution P4. The four plots correspond to t = Q^T/A^T/2^2>T/A^ where T is the period of P4. The 
plots are in clockwise order, beginning at the top left. The velocity fields were shifted in such a 
way that a plot at t = T would coincide exactly with the plot for t = 0. The isolines are thin if 
the streamwise velocity is positive, and thick if it is negative. Each of the four plots has 24 isolines 
equispaced between minimum values which are —0.5, —1.4, —2.3, —0.5 (in wall units) and maximum 
values which are 9.4, 9.4, 7.5, 8.3. 
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Figure 7: The plot above shows the variation of the root mean square value of the fluctuating 
part of streamwise velocity, which is denote by u'~^, in the spanwise direction. The rms value is 
computed at y^ ~ 10 and x = ttAx- 



coordinate but varies much more strongly in the spanwise direction. In contour plots such as those 
in Figure m streakiness shows up as isolines that are nearly parallel to the x axis. It is clear from 
that figure that streaks at i = T/2 break-up at t = 3r/4 and then re-form. The plot at t = 
corresponds to the initial velocity field of P4. 

It is difficult to measure the velocity field as a whole in experiments, and therefore experimental 
visualizations of s treaks rely on hyd r ogen bubbles introduced i nto the flow by platinum wires and 
other techniques (JKline et al.l . 119671: ISmith and Metzleii . Il983l ). Sometimes streaks are detected 
using pointwise measurements (JKlebanoff et al.l . Il962l ). Figure [7] shows the dependence of the 
fluctuations of the streamwise velocity on the spanwise direction. The rms velocity n'"*" shown in 
the top-right plot of Figure [5] was obtained by averaging over a period and by averaging spatially 
in the spanwise and streamwise directions. But in Figure [7l only the time average was used for 
computing the rms velocity. 

The possibility that the break-up of streaks o bserved experiinental ly could be an artifact of flow 
visualization techniques has been considered by I Jimenez et al.l (120051 ). It has been suggested that 
the advection of permanent objects could be partly responsible for experimental observations. Our 
computation of relative periodic solutions (P2 through P5) shows that temporal periodicity and 
the advection of permanent objects are not mutually exclusive possibilities. While it is generally 
implicitly assumed that the advection of coherent structures could only be in the streamwise direc- 
tion, our computation of P2 and P4 shows that the advection could be in the spanwise direction as 
well. This could be significant, as will be seen shortly. 

The spanwise variation in the strength of the fiuctuations of the streamwise velocity shown in 



Figur e [71 for Pi through Pg is reminiscent of experimental data, such as Figure 2 of (JKlebanoff et al 



19621 ) for instance. This spanwise variation is not very pronounced for Pi, which is not a bursting 
solution, but is very pronounced for P2 and P4. Those are the only two solutions in Table [1] which 
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Figure 8: The two plots above show slices through the streamwise velocity field of Pi at t = and 
t = T/2, respectively. The slices were taken at y"*" ~ 9. As explained in the text, the plots show 
that Pi does not have the shift-reflection symmetry. Isolines of u~^ are drawn at 24 equispaced 
values between a maximum of 11.6 and a minimum of 1.6. The maximum occurs in the wide gap 
between isolines in the lower part of either plot, and the minimum occurs in the gap in the upper 
part. 

have a spanwise drift. We are led to conclude that spanwise advection of coherent structures could 
be a significant source of the observed spanwise variation of n'+. 



4-4 Discrete symmetries of plane Couette flow 



The Navier-Stokes equation for plane Couette flow has two discrete symmetries (jWaleffd . l2003l ) 
The shift-reflection transformation of the velocity field is given by 





and the shift-rotation transformation of the velocity field is given by 



-x + TTAx,-y,z -hvrA^ 



Plane Couette flow is unchanged under both these transformations. Thus if a single velocity field 
along a trajectory of plane Couette flow satisfies either symmetry, all points along the trajectory 
must have the same symmetries. However, velocity fields that lie on the stable and unstable 
manifolds of symmetric periodic or relative periodic solutions need not be symmetric. 

Figure [8] shows that Pi does not have the shift-reflection symmetry. Both the plots would look 
very different if they were flipped upside-down and shifted forward by half the width of the plot 
in the x direction. We have verified that Pi does not have the shift-rotation symmetry either. A 
close inspection of Figure [8] reveals that the second plot can be obtained by shifting the first plot 
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Figure 9: The plot above shows the projection of the unstable manifold of Pi to the energy dissi- 
pation vs. energy input plane. 



in the x direction. In fact, Pi is a relative periodic solution. During half the period listed in Table 
[U the initial velocity field of Pi shifts by Sx/ttAx = 0.5 and Sz = 0. 

The averaged velocity field of a long turbulent trajectory of plane Couette flow satisfies the 
discrete symmetries to a very rough approximation. However, there is no reason to believe that 
such an average has any dynamical significance. The average could just be a transient state and it 
may not even lie on the asymptotic set of the flow. 

The six solutions listed in Table [1] satisfy neither the shift-reflection symmetry nor the shift- 
rotation symmetry. This implies that the shift-reflection and shift-rotation transformations (which 
commute with each other) can be applied to each solution in Table [1] to obtain 18 more solutions. 

4-5 Unstable manifolds of the solutions Pi 



In their study of the bursting phenomenon in turbulent Poiseuille flow, lltano and Tohl (120011 ) 
found that the break-up and re-formation of coherent structures is well approximated by the un- 
stable manifold of a traveling wave sol ution. Among t he so lutions we have computed. Pi most 
resembles the traveling wave solution of (lltano and Tohl . l200ll ) . Motivated by that resemblence, we 
asked if the solutions P2 through Pg could be related to the unstable manifold of Pi . 

We perturbed Pi along the single unstable direction that corresponds to that solution. A 
perturbation in one sense begins to head towards the laminar solution right away and falls into 
the laminar solution. A perturbation with the opposite sign heads into the turbulent region, as 
shown in Figure [9| but relaminarizes eventually. In contrast, we found that perturbations along 
the unstable manifolds of the bursting periodic solutions P2 through Pg do not relaminarize. Thus 
there appears to be no relationship between the unstable manifold of Pi and the other solutions. 

In addition, unlike Pi the bursting solutions have more than one unstable direction. For in- 
stance, the Arnoldi iteration shows that P2 has at least 11 unstable directions. These facts suggest 
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that the unstable manifold of Pi cannot by itself explain the bursting phenomenon in the plane 
Couette scenario considered in this paper. 

5 Discussion 

The approach to the bursting phenomenon, and to turbulence in fluid flows more generally, that 
this paper advocates is quite simple. As the incompressible Navier-Stokes equation is an excellent 
physical model, our approach is simply to compute solutions of that differential equation. This ap- 
proach should not be surprising in itself because the main reason for deriving differential equations 
is to understand their solutions. However, well-resolved spatial discretizations of turbulent phe- 
nomena require more than a hundred thousand degrees of freedom and computing solutions with 
so many degrees of freedom is nontrivial. We have shown that problems associated with largeness 
of the number of degrees of freedom can be overcome. Indeed, our computations of periodic and 
relative periodic solutions used as many as twenty times the number of degrees of freedom in any 
earlier computation of periodic solutions. 

Having decided to compute solutions, one must decide what type of solutions to look for within 
turbulent flows. Here the basic principles of dynamics are of help. The recurrent nature of the 
bursting phenomenon has been noted by experimentalists, as we pointed out earlier. Recurrence, 
which is the tendency of certain dynamical systems to revisit points close to their initial state 
in phase space after extended excursions, has long been a central and unifying idea in dynamics. 
The relevance of recurrence to long term dynamics is cap tured in a most general way by the 
Poincare recurrence theorem ( Katok and Hasselblatd . Il995l ). In its early days, this theorem gave 



rise to troubling questions about the foundations of statistical mechanics. The number of degrees 
of freedom needed for resolving low Reynolds number turbulence is quite large, yet far smaller than 
the number of degrees of freedom typical of statistical mechanics. Another significant difference is 
that low Reynolds number turbulence is governed by equations that are not Hamiltonian and in 
which viscosity damps high wavenumbers. Therefore there is reason to think that a point of view 
based on recurrences will be useful for understanding turbulence. 

A combination of local instability, which is a well-known feature of turbulent flows, and bound- 
edness in phase space naturally suggests the existence of periodic motions. The closing lemma 
is a princi ple that suggests that dynarn ics can be understood quite generally in terms of periodic 
solutions ( Katok and Hasselblattl . Il995l ) . It has been proved in a few restricted settings and serves 



as a beacon. While considering the closing lemma, it is worth noting again that the right notion 
of recurrence depends upon invariance properties of the underlying differential equation. If the 
differential equation is unchanged by a continuous group of transformations, then it is appropriate 
to look for relative periodic solutions. 

In som e well-unders tood settings, it is possible to prove the existence of infinitely many periodic 
solutions ( Moseii . Il973l ). Even though turbulent phenomena are not known to fall under any of 



these settings, there is reason to think that there are infinitely many periodic and relative periodic 
solutions embedded within turbulent flows. The lack of complete theoretical results should not be an 
impediment to computation, since the major purpose of computation is to render tractable problems 
that are beyond the reach of theory. An advantage of computing many of these solutions is that 
they can help us understand the dynamics as a whole in terms of accurately computed periodic 
solutions. In addition, these periodic s olutions could serve as a basis to understand turbulent 
statistics using the periodic orbit theory ( Cvitanovic et al.l . l2005l ). 
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It is natural to expect different solutions of the same differential equation to bear a relationship 
to one another. In the case of linear differential equations, the principle of linear superposition 
gives this relationship. For chaotic nonlinear systems, the relationship between solutions is much 
more complicated and intrig uing. I i i seve ral instances, the precise nature of the relationship is 



given by symbolic dynamics (IMoseii. Il973l). Algor i thms for computing nonlinear systems can be 



based upon symbolic dynamics ( Viswanathl . l2003l . 1200J). Those algorithms made it possible to 



compute the fractal structure of the well- known Lorenz attractor. Although the fractal structure 
of the Lorenz attractor was deduced by iLorend (jl963l ). its computation became possible only 



with algorithms based on symbolic dynamics. The Lorenz equations were derived to illuminate 
the essentially deterministic nature of turbulence and they have been completely successful in that 
respect. Computations of the strange attractor of the Lorenz equations can likewise serve as models 
for computations of turbulence. 

The Lorenz computations based on symbolic dynamics illustrate the advantanges of periodic 
solutions over steady solutions in understanding the dynamics as a whole. It is possible to accurately 
compute periodic solutions that are as close as machine precision permits to ra ndom points on th e 
Lorenz attractor, and thus obtain a very good understanding of the dynamics ( Viswanathl . |2003| ). 



Such a precise understanding cannot be obtained using steady solutions alone. To understand 
important aspects of turbulent boundary layers such as the turbulent energy production in the 
buffer layer, it is likewise necessary to go beyond steady solutions and traveling waves and compute 
periodic and relative periodic solutions. 

It would be premature to suggest that the various solutions embedded within turbulent flows 
can be described using symbolic dynamics. Yet, it would be very surprising if these solutions did 
not bear a relationship to one another. Discovering such a relationship would be a major advance 
in our understanding of turbulent flows. 

This paper has asserted the existence of six solutions of turbulent plane Couette flow. Error 
estimates for these solutions were supplied in Tabled! We argued in Section 3 that these solutions, 
along with their error estimates, can be verified using a good code for direct numerical simulation 
of channel flows. Such quantitative reproducibility is a step forward in the area of turbulence 
computation. 
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